{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Blind kildeseparasjon med uavhengig komponentanalyse\n",
    "\n",
    "Prosjekt 1 i TMA4320 våren 2019 - av Thorvald Ballestad, Jonas Bueie og Herman Sletmoen (gruppe 3)\n",
    "\n",
    "*Blind kildeseparasjon* (engelsk: *blind source separation*) er separasjon av en mengde kildesignaler fra observerte blandinger av disse, uten informasjon om kildesignalene og blandingsprosessen.\n",
    "\n",
    "*Cocktailparty-problemet* er et konkret eksempel på blind kildeseparasjon.\n",
    "Betrakt et område med $n$ lydkilder som ved tiden $t$ spiller av lydsignaler med amplituder $s_1(t), \\ldots, s_n(t)$, samtidig som $n$ mikrofoner tar opp lydsignaler med amplituder $x_1(t), \\ldots, x_n(t)$.\n",
    "Mikrofonene registrerer totale signaler $x_i(t) = a_{i1} s_1(t) + \\ldots + a_{in} s_n(t)$ som lineære kombinasjoner av kildesignalene, der $a_{ij}$ er en antatt tidsuavhengig faktor som bestemmer bidraget fra kildesignalet $s_j(t)$ til totalsignalet $x_i(t)$.\n",
    "Faktorene $a_{ij}$ avhenger gjerne av geometri, og vil i et perfekt fysisk scenario være omvendt proporsjonale med kvadratet av avstanden $r_{ij}$ mellom lydkilden $s_j(t)$ og mikrofonen $x_i(t)$, som kjent fra klassisk fysikk.\n",
    "\n",
    "Med kolonnevektorene $\\boldsymbol{x} = [x_1(t), \\ldots, x_n(t)]^T$ og $\\boldsymbol{s} = [s_1(t), \\ldots, s_n(t)]^T$ og *miksematrisen* $A = [a_{ij}]$ av størrelse $n \\times n$, kan vi mer konsist skrive\n",
    "\n",
    "$$\\boldsymbol{x} = A \\boldsymbol{s}$$\n",
    "\n",
    "I en typisk situasjon vil mikrofonene sample totalsignalene $x_i(t_j)$ ved $N$ diskrete tidspunkter $t_1, \\ldots, t_N$, og vi ønsker å bestemme kildesignalene $s_i(t_j)$ ved de samme tidspunktene som bygger opp totalsignalet.\n",
    "Vektorligningen over vil da gjelde ved ethvert tidspunkt $t_j$.\n",
    "Med totalsignalsmatrisen $X = [x_{ij}] = [x_i(t_j)]$ og kildsignalsmatrisen $S = [s_{ij}] = [s_i(t_j)]$, begge av størrelse $n \\times N$, kan vi dermed oppsummere med matriseligningen\n",
    "\n",
    "$$X = A S$$\n",
    "\n",
    "I cocktailparty-problemet er kun de observerte totalsignalene $X$ kjent, og vi søker en invers matrise $W$, gjerne kalt *demiksematrisen*, til miksematrisen $A$ slik at $W X$ gir oss de opprinnelige kildesignalene.\n",
    "Oppgaven virker tilsynelatende umulig, da vi hverken har kunnskap om $S$ eller $A$.\n",
    "I praksis, derimot, viser det seg at vi kan oppnå overraskende gode estimater for $S$ gjennom en enkel antagelse om den statistiske distribusjonen av kildesignalene.\n",
    "\n",
    "Vi vil her løse cocktail-party problemet med *uavhengig komponentanalyse* (engelsk: *independent component analysis (ICA)*).\n",
    "Kjerneantagelsen i ICA er at de opprinnelige signalkomponentene $s_i(t)$ er **statistisk uavhengige**.\n",
    "For å grundig forstå hvordan denne antagelsen legger grunnlaget for ICA-metoden, må vi se den i sammenheng med sentralgrenseteoremet fra statistikken.\n",
    "\n",
    "La oss definere de stokastiske variablene $X_i$ og $S_i$ med fordelinger som representerer et utvalg av totalsignalene $x_i(t_j)$ og kildesignalene $s_i(t_j)$, henholdsvis.\n",
    "Vi definerer så vektorene $\\boldsymbol{X} = [X_1, \\ldots, X_n]^T$ og $\\boldsymbol{S} = [S_1, \\ldots, S_n]^T$, slik at sammenhengen mellom distribusjonene blir $\\boldsymbol{X} = A \\boldsymbol{S}$.\n",
    "Vi ønsker da å finne matrisen $W$ slik at $W \\boldsymbol{X}$ gir oss tilbake distribusjonene av kildesignalene.\n",
    "Hvilken egenskap må en rad $\\boldsymbol{w} = [w_1, \\ldots, w_n]$ i $W$ ha, gitt at kildesignalene $S_i$ er statistisk uavhengige?\n",
    "\n",
    "La oss nå betrakte den stokastiske variabelen $Y = \\sum_i w_i \\cdot X_i = \\boldsymbol{w} \\boldsymbol{X}$.\n",
    "Siden $\\boldsymbol{w}$ er en rad i $W$, vil $Y$ representere et av kildesignalene $S_i$.\n",
    "Ved å definere $\\boldsymbol{z} = [z_1, \\ldots, z_n] = \\boldsymbol{w} A$, kan vi skrive $Y = \\boldsymbol{w} A \\boldsymbol{S} = \\boldsymbol{z} \\boldsymbol{S} = \\sum_i z_i \\cdot S_i$.\n",
    "Altså er $Y$ en lineær kombinasjon av de uavhengige stokastiske variablene $S_1, \\ldots, S_n$.\n",
    "Fra sentralgrenseteoremet vet vi at en lineær kombinasjon av to eller flere uavhengige stokastiske variabler sannsynligvis er mer gaussisk fordelt enn de uavhengige variablene.\n",
    "Dermed er $Y$ minst like gaussisk fordelt som $S_i$-ene, og den blir minst gaussisk fordelt når kun ett av elementene i $\\boldsymbol{z}$ er forskjellig fra $0$.\n",
    "For at $Y$ skal representere fordelingen til et av kildesignalene $S_i$, må vi derfor velge $\\boldsymbol{w}$ slik at $\\boldsymbol{w} \\boldsymbol{X}$ er minst mulig gaussisk fordelt!\n",
    "\n",
    "I denne notatboken vil vi implementere *FastICA*-algoritmen, utarbeidet av Aapo Hyvärinen ved Helsinki University of Technology for å løse cocktailparty-problemet.\n",
    "\n",
    "FastICA-algoritmen kan oppsummeres **svært grovt** i følgende steg:\n",
    "1. Preprossesér dataene $X$ ved å sentrere dem.\n",
    "2. Preprossesér dataene $X$ ved å white dem.\n",
    "3. Gjett en verdi $W_0$ for demiksematrisen.\n",
    "4. Finn en bedre approksimasjon $W_k$ for demiksematrisen for $k = 1, 2, \\ldots$ som gjør radene i $W_k X$ mer gaussisk fordelt enn i $W_{k-1} X$, inntil $W_k = W_{k-1} = W$.\n",
    "5. Beregn de antatte kildesignalene gjennom transformasjonen $W X$ av totalsignalene.\n",
    "\n",
    "Vi vil ha en anvendt tilnærming der vi motiverer og forklarer hensikten bak de ulike stegene i algoritmen mens vi implementerer den del for del, som beskrevet i et dokument av Brynjulf Owren.<sup>[1]</sup>\n",
    "For detaljerte bevis for \"korrektheten\" til algoritmen, refererer vi til materiale produsert av Hyvärinen selv.<sup>[2]</sup>\n",
    "Til slutt vil vi teste algoritmen på noen lydklipp."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I denne notatboken benytter vi noen ulike Python-biblioteker:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np # for effective numerical operations\n",
    "import wav_file_loader # for loading data from .wav audio files\n",
    "import IPython.display as ipd # for interactive audio playback\n",
    "import time # for benchmarking"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`wav_file_loader` vil utleveres sammen med notatboken, mens `numpy` og `IPython` anses som velkjente og lett tilgjengelige biblioteker."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprosessering\n",
    "\n",
    "I FastICA-algoritmen vil dataene $X$ først sendes gjennom et par preprosseseringssteg.\n",
    "Hensikten med preprosseseringen er å gi dataene egenskaper som forenkler senere steg i algoritmen.\n",
    "\n",
    "Først vil amplitudene i hvert lydopptak forskyves gjennom\n",
    "\n",
    "$$\\boldsymbol{x_i} \\gets \\boldsymbol{x_i} - E[\\boldsymbol{x_i}]$$\n",
    "\n",
    "der $\\boldsymbol{x_i}$ er rad $i$ i $X$ og $E[\\boldsymbol{x_i}]$ er forventningsverdien til elementene i raden, slik at dataene i radene er sentrert om sine middelverdier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def center_row_means(M):\n",
    "    \"\"\"Returns a new matrix in which the elements of each row in\n",
    "    the matrix M are shifted so the mean of each row is zero.\n",
    "    \"\"\"\n",
    "    n_rows = len(M)\n",
    "    means = np.mean(M, axis=1).reshape((n_rows, 1))\n",
    "    M -= means # subtract the mean of each row from all its elements\n",
    "    return M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Deretter vil dataene sendes gjennom \"whitening\"-transformasjonen\n",
    "\n",
    "$$ X \\gets C^{-1/2} X$$\n",
    "\n",
    "der $C = \\frac{1}{N - 1} X X^T$ er en såkalt kovariansematrise til dataene $X$.\n",
    "Med $C^{-1/2}$ mener vi en matrise som oppfyller $C^{-1/2} C^{-1/2} = C^{-1}$.\n",
    "Vi finner denne såkalte kvadratroten til inversmatrisen ved hjelp av egenverdidekomponeringen $C^{-1/2} = E D^{-1/2} E^{-1}$, der $E$ er egenverdimatrisen til $C$ og $D^{-1/2}$ er en diagonal matrise med korresponderende inverse kvadratrøtter av egenverdiene til $C$ langs diagonalen.\n",
    "Konstruksjonen av $C$ som en symmetrisk matrise garanterer eksistensen av reelle og positive egenverdier og gir i tillegg forenklingen $E^{-1} = E^T$.\n",
    "\n",
    "Whitening-transformasjonen endrer dataene $X$ på en slik måte at ligningen $X = A S$ fortsatt skal løses under de samme forutsetningene som før, men mot en **ortogonal** matrise $A$, noe som gir ytterligere forenklinger senere i algoritmen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def symmetric_matrix_inverse_sqrt(M):\n",
    "    \"\"\"Returns a new matrix that is the square root of the inverse of the matrix M.\n",
    "    The matrix M is assumed to be invertible and symmetric.\n",
    "    \"\"\"\n",
    "    eigvals, eigvecs = np.linalg.eig(M)\n",
    "    D = np.diag(1 / np.sqrt(eigvals))\n",
    "    M_inv_sqrt = eigvecs @ D @ eigvecs.T\n",
    "    return M_inv_sqrt\n",
    "\n",
    "def whiten_matrix(X):\n",
    "    \"\"\"Returns a new whitened matrix from the matrix X.\"\"\"\n",
    "    n_cols = X.shape[1]\n",
    "    C = (1 / (n_cols - 1)) * X @ X.T\n",
    "    C_inv_sqrt  = symmetric_matrix_inverse_sqrt(C)\n",
    "    X = C_inv_sqrt @ X \n",
    "    return X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iterasjon for demiksematrisen\n",
    "\n",
    "Kjernen i FastICA-algoritmen er iterasjonen over matriser $W_k$ som maksimerer ikke-gaussiskheten til radene i $W_k X$.\n",
    "Iterasjonen bestemmer en bedre approksimasjon $W_k$ til $W$ basert på en tidligere approksimasjon $W_{k-1}$, slik at radene i $W_k X$ er mer gaussisk fordelt enn i $W_{k-1} X$.\n",
    "Vi understreker at dataene $X$ på dette stadiet i algoritmen er preprossesert i henhold til stegene over.\n",
    "\n",
    "Først måles gaussiskheten til dataene i radene til $S_{k-1} = W_{k-1} X$ ved å la en velvalgt reell funksjon $g(u)$ og dens deriverte $g'(u)$ operere elementvis på $S_{k-1}$.\n",
    "Resultatene lagres i matrisene $G$ og $G'$.\n",
    "Vi åpner her for bruk av både kurtosefunksjonen $g_1(u) = 4 u^3$ og negentropifunksjonen $g_2(u) = u e^{-u^2/2}$.\n",
    "Matrisene $G$ og $G'$ inneholder dermed informasjon om den \"nåværende\" gaussiskheten til dataene og \"endringen\" i gaussiskheten ved variasjon av elementer i matrisen $W_k$.\n",
    "Basert på denne informasjonen, beregnes først en ny matrise\n",
    "\n",
    "$$W_k^+ = \\frac{1}{N} G X^T - \\text{diag}(E[G']) W_{k-1}$$\n",
    "\n",
    "der $\\text{diag}(E[G'])$ er en diagonal $n \\times n$ matrise med middelverdiene til radene i $G'$ langs diagonalen.\n",
    "Så normaliseres radene i $W_k^+$ slik at den euklidske normen av hver rad er $1$.\n",
    "\n",
    "Til slutt beregnes $W_k$ som den nærmeste ortogonale matrisen til $W_k^+$ gjennom transformasjonen\n",
    "$$W_k = (W_k^+ (W_k^+)^T)^{-1/2} W_k^+$$\n",
    "\n",
    "Kvadratroten av en invers matrise beregnes her på samme måte og under samme forutsetninger som før.\n",
    "\n",
    "Vi har her delt opp denne prosedyren i et optimaliseringssteg og et ortogonaliseringssteg, også kjent som et dekorreleringssteg."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def g_kurtosis(u):\n",
    "    \"\"\"Kurtosis\"\"\"\n",
    "    return 4 * u**3\n",
    "\n",
    "def dg_kurtosis(u):\n",
    "    \"\"\"Kurtosis derivative\"\"\"\n",
    "    return 12 * u**2\n",
    "\n",
    "def g_negentropy(u):\n",
    "    \"\"\"Negentropy\"\"\"\n",
    "    return u * np.exp( -u**2 / 2)\n",
    "    \n",
    "def dg_negentropy(u):\n",
    "    \"\"\"Negentropy derivative\"\"\"\n",
    "    return (1 - u**2) * np.exp( -u**2 / 2)    \n",
    "\n",
    "def normalize_row_norms(M):\n",
    "    \"\"\"Returns a new matrix in which the rows of the matrix M \n",
    "    are normalized so their euclidean norm is 1.\n",
    "    \"\"\"\n",
    "    n_rows = len(M)\n",
    "    row_norms = np.linalg.norm(M, axis=1).reshape((n_rows, 1))\n",
    "    M /= row_norms # divide each row by its norm\n",
    "    return M\n",
    "\n",
    "def optimize_demixing_matrix(W, X, g, dg):\n",
    "    \"\"\"Returns a new matrix W' so that the rows of W' * X are distributed more non-gaussianly\n",
    "    than the rows of W * X.\n",
    "    \"Gaussianity\" is measured using the scalar function g and its derivative, dg.\n",
    "    The rows of the matrix X are assumed to have zero mean, and X is assumed to be whitened.\n",
    "    \"\"\"\n",
    "    S = W @ X\n",
    "    G = g(S)\n",
    "    dG = dg(S)\n",
    "    N = X.shape[1]\n",
    "    W = (1 / N) * G @ X.T - np.diag(np.mean(dG, axis=1)) @ W\n",
    "    W = normalize_row_norms(W)\n",
    "    return W\n",
    "\n",
    "def orthogonalize_matrix(W):\n",
    "    \"\"\"Returns a new matrix that is the nearest orthogonal matrix to the matrix W.\"\"\"\n",
    "    return symmetric_matrix_inverse_sqrt(W @ W.T) @ W\n",
    "\n",
    "def update_demixing_matrix(W, X, g, dg):\n",
    "    \"\"\"Returns a new matrix that is the next iteration of the demixing matrix in the FastICA algorithm,\n",
    "    given the previous demixing matrix W and the data matrix X.\n",
    "    The rows of the matrix X are assumed to have zero mean, and X is assumed to be whitened.\n",
    "    \"\"\"\n",
    "    W = optimize_demixing_matrix(W, X, g, dg)\n",
    "    W = orthogonalize_matrix(W)\n",
    "    return W"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FastICA-algoritmen"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nå som vi har funksjonalitet for å finne stadig bedre approksimasjoner til demiksematrisen $W$, er vi endelig klare til å sette sammen alle delene til FastICA-algoritmen som beskrevet over!\n",
    "\n",
    "Algoritmen itererer over stadig bedre approksimasjoner $W_k$ for $W$ inntil to etterfølgende matriser $W_k$ og $W_{k-1}$ har konvergert innenfor en gitt toleranse eller iterasjonen er gjennomført et gitt maksimalt antall ganger.\n",
    "Konvergens måles her med det maksimale avviket på diagonalen mellom identitetsmatrisen $I$ og $W_k W_{k-1}^T$.\n",
    "Dette er mer effektivt i denne sammenhengen enn å, for eksempel, finne det maksimale elementvise avviket mellom $W_k$ og $W_{k-1}$, da algoritmen garanterer at $W_k$ og $W_{k-1}$ alltid er ortogonale, slik at $W_k W_k^T = W_{k-1} W_{k-1}^T = I$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def orthonormal_matrices_are_equal(A, B, tol):\n",
    "    \"\"\"Returns whether the matrices A and B are equal up to the tolerance tol. \n",
    "    A and B are assumed to be orthonormal matrices.\n",
    "    The tolerance is compared with the maximum deviance between\n",
    "    the identity matrix and A * B^T along their diagonals.\n",
    "    \"\"\"\n",
    "    row_dot_products = np.sum(A * B, axis=1)\n",
    "    max_dev = np.abs(np.max(1 - np.abs(row_dot_products)))\n",
    "    return max_dev < tol\n",
    "\n",
    "def fast_ICA(X, tol=1e-10, max_iters=100, g=g_kurtosis, dg=dg_kurtosis):\n",
    "    \"\"\"Returns new matrices S and W containing separated data components and the \n",
    "    demixing matrix W using the FastICA algorithm, given data from the matrix X,\n",
    "    and the number of iterations, i, the algorithm uses to approximate W.\n",
    "    The algorithm terminates after reaching max_iters iterations or when \n",
    "    two successive approximations of W are equal within the tolerance tol.\n",
    "    \"\"\"\n",
    "    n = len(X)\n",
    "    X = center_row_means(X)\n",
    "    X = whiten_matrix(X)\n",
    "    \n",
    "    W = np.identity(n) # initial guess for W\n",
    "    for i in range(0, max_iters):\n",
    "        W0 = W\n",
    "        W = update_demixing_matrix(W, X, g, dg) \n",
    "        if orthonormal_matrices_are_equal(W, W0, tol):\n",
    "            break\n",
    "\n",
    "    if i == max_iters - 1:\n",
    "        print(\"warning: FastICA returned before convergence\")\n",
    "    S = W @ X\n",
    "    return S, W, i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Testing av FastICA-algoritmen\n",
    "\n",
    "I denne seksjonen anvender vi FastICA-algoritmen på flere blandede lydopptak og sammenligner de separerte lydsignalene med de blandede opptakene og eventuelle kjente, opprinnelige kildesignaler.\n",
    "\n",
    "For å forenkle denne testeprosessen, har vi utviklet et enkelt rammeverk som håndterer innlesing og avspilling av lyddata samt operasjoner på disse.\n",
    "Rammeverket vil normalisere både blandede og separerte lydsignaler, slik at de oppfattes til å ha lik maksimal styrke.\n",
    "Det har også funksjonalitet for å blande sammen lydsignaler på en tilfeldig måte.\n",
    "\n",
    "Detaljene i rammeverket er på ingen måte nødvendig å undersøke for andre enn spesielt interesserte, da dets eneste hensikt er å forenkle og tydeliggjøre testingen av FastICA-algoritmen på forskjellige lydfiler og på ulike måter.\n",
    "Hvordan testingen foregår vil gå klart frem av bruken av rammeverket.\n",
    "Vi understreker allikevel for eventuelle skeptiske lesere at rammeverket alltid separerer blandede lydklipp med FastICA-algoritmen, enten de opprinnelige separerte lydklippene er kjente eller ikke."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def normalize_audio(data):\n",
    "    \"\"\"Returns a new matrix in which the rows of the matrix M \n",
    "    are normalized so their maximum absolute values are 1.\n",
    "    \"\"\"\n",
    "    n_tracks = len(data)\n",
    "    row_maximums = np.amax(np.abs(data), axis=1).reshape((n_tracks, 1))\n",
    "    data /= row_maximums # divide each row by its maximum absolute value\n",
    "    return data\n",
    "\n",
    "def random_mixing_matrix(signals, observations):\n",
    "    \"\"\" Returns a new random mixing matrix.\n",
    "    Each element is a small positive number, not too close to 0.\n",
    "    \"\"\"\n",
    "    A = 0.25 + np.random.rand(observations, signals)\n",
    "    return normalize_audio(A)\n",
    "\n",
    "class Tester:\n",
    "    \"\"\"Object used to demonstrate and test the FastICA algorithm \n",
    "    on the .wav audio files located at given file paths.\n",
    "    \"\"\"\n",
    "    def __init__(self, paths):\n",
    "        data_observed, srate = wav_file_loader.read_wavefiles(paths)\n",
    "        self.data_observed = normalize_audio(data_observed)\n",
    "        self.data_separated = None\n",
    "        self.srate = srate\n",
    "        self.paths = paths\n",
    "        self.set_tol(1e-9)\n",
    "        self.use_kurtosis()\n",
    "\n",
    "    def mix_observed(self):\n",
    "        print(\"Mixing observed data\")\n",
    "        mixmat = random_mixing_matrix(len(self.paths), len(self.paths))\n",
    "        self.data_observed = mixmat @ self.data_observed\n",
    "        self.data_observed = normalize_audio(self.data_observed)\n",
    "        \n",
    "    def set_tol(self, tol):\n",
    "        self.tol = tol\n",
    "    \n",
    "    def set_g(self, g):\n",
    "        self.g = g[0]\n",
    "        self.dg = g[1]\n",
    "    \n",
    "    def use_kurtosis(self):\n",
    "        \"\"\"Assumes g_kurtosis and dg_kurtosis defined.\"\"\"\n",
    "        self.set_g((g_kurtosis, dg_kurtosis))\n",
    "        \n",
    "    def use_negentropy(self):\n",
    "        \"\"\"Asssumes g_negentropy and dg_negentropy defined.\"\"\"\n",
    "        self.set_g((g_negentropy, dg_negentropy))\n",
    "        \n",
    "    def play_observed(self):\n",
    "        print(\"Playing observed data:\")\n",
    "        self._play_data(self.data_observed)\n",
    "        \n",
    "    def play_separated(self, sample_percent_of_data=1):\n",
    "        if sample_percent_of_data == 1:\n",
    "            self._separate_data()\n",
    "        else:\n",
    "            # create sampler that selects a random sample of data\n",
    "            sampler = lambda data_list: data_list[:, np.random.randint(data_list.shape[1],size = int(sample_percent_of_data*data_list.shape[1]))]\n",
    "            self._separate_data(sampler)\n",
    "            \n",
    "        print(\"Playing separated data:\")\n",
    "        self._play_data(self.data_separated)\n",
    "        \n",
    "    def get_duration(self):\n",
    "        return self.data_observed.shape[1] / self.srate\n",
    "        \n",
    "    def _play_data(self, data_list):\n",
    "        for data in data_list:\n",
    "            ipd.display(ipd.Audio(data, rate=self.srate))\n",
    "    \n",
    "    def _separate_data(self, sampler = lambda x: x):\n",
    "        start = time.time()\n",
    "        \n",
    "        data_observed_sample = sampler(self.data_observed)\n",
    "        percent_of_data = (data_observed_sample.shape[1] / self.data_observed.shape[1]) * 100\n",
    "        \n",
    "        print(\"Started separating data\")\n",
    "        print(\"  Number of tracks:    %d\" % self.data_observed.shape[0])\n",
    "        print(\"  Samples per track:   %d\" % self.data_observed.shape[1])\n",
    "        print(\"  Sample rate:         %d Hz\" % self.srate)\n",
    "        print(\"  Duration per track:  %.1f s\" % self.get_duration())\n",
    "        print(\"  Gaussianity measure: %s\" % self.g.__doc__.split()[0].lower())\n",
    "        print(\"  Tolerance:           %e\" % self.tol)\n",
    "        print(\"  Sample size:         %.2f%%\" % percent_of_data)\n",
    "        \n",
    "        # do transformation of data ourselves\n",
    "        _, W, iters = fast_ICA(data_observed_sample, tol=self.tol, g=self.g, dg=self.dg)  \n",
    "        self.data_separated = whiten_matrix(center_row_means(self.data_observed))\n",
    "        self.data_separated = W @ self.data_separated\n",
    "        self.data_separated = normalize_audio(self.data_separated)\n",
    "        \n",
    "        end = time.time()\n",
    "        \n",
    "        print(\"Finished separating data\")\n",
    "        print(\"  Iterations:          %d\" % iters)\n",
    "        print(\"  Time used:           %.4f s\" % (end - start))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vi tester først FastICA-algoritmen på $n = 3$ utdelte blandede lydopptak av like mange lydkilder.\n",
    "Testen gjøres både med kurtosevarianten og negentropivarianten av funksjonene som måler gaussiskhet.\n",
    "\n",
    "Vi hører at lydsignalene separeres på en svært god måte, både ved bruk av kurtose og negentropi i målingen av gaussiskheten til data underveis i algoritmen.\n",
    "Alle de separerte lydsporene er helt rene og meningsfulle lydspor og ingen av dem bærer preg av å være blandet med noen av de andre separerte lydsporene.\n",
    "Lydsporene separeres like godt uavhengig av om de inneholder tale eller musikk."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "paths = [\"audio/\" + f + \".wav\" for f in (\"mix_1\", \"mix_2\", \"mix_3\")]\n",
    "test = Tester(paths)\n",
    "\n",
    "test.play_observed()\n",
    "\n",
    "test.use_kurtosis()\n",
    "test.play_separated()\n",
    "\n",
    "test.use_negentropy()\n",
    "test.play_separated()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Vi illustrerer nå hvordan langvarige blandede lydklipp kan separeres svært effektivt ved å kun bruke et tilfeldig **utvalg** av totalsignalene i beregningen av demiksematrisen $W$, for å så benytte denne matrisen til å separere **hele** det blandede lydklippet.\n",
    "Den naive metoden med å benytte hele lydklippet også i beregningen av $W$ fungerer like godt, men ikke merkbart bedre, og er betydelig tregere og krever mer minne.\n",
    "Leseren kan selv sammenligne tidsbruken med den naive metoden på de samme lydklippene ved å fjerne kommentaren av linjen under som gjennomfører den naive metoden.\n",
    "\n",
    "Intuitivt er det tenkelig at miksematrisen $A$, og dermed demiksematrisen $W$, avhenger av romlige og andre tidsuavhengige faktorer.\n",
    "Dermed er det interessant, men ikke overraskende, at kun et utvalg av totalsignalene er nødvendig i bestemmelsen av $W$, da prosedyren bygger på statistiske prinsipper om distribusjoner av data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "paths = [\"audio/\" + f + \".wav\" for f in (\"nwa1\", \"nwa2\", \"nwa3\", \"jimcarreymb\")]\n",
    "test_mix = Tester(paths)\n",
    "\n",
    "test_mix.play_observed()\n",
    "test_mix.mix_observed()\n",
    "test_mix.play_observed()\n",
    "#test_mix.play_separated() # commented by default\n",
    "test_mix.play_separated(sample_percent_of_data=0.02)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Til slutt forsøker vi å benytte FastICA-algoritmen på flere, kortere lydklipp.\n",
    "Vi benytter her 8 ulike lydklipp.\n",
    "Her får vi best resultater om vi benytter negentropi og bruker hele lydklippene for å beregne demiksematrisen.\n",
    "\n",
    "Vi oppfordrer leseren til å forsøke å benytte kurtose og å separere dataene kun med halvparten av datapunktene ved å fjerne kommentarene på en av de kommenterte linjene under eller begge.\n",
    "I disse tilfellene overlapper flere av de separerte lydsporene på vår maskin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "paths = [\"audio/\" + f + \"_short.wav\" for f in (\"nwa1\", \"nwa2\", \"nwa3\", \"nwa4\", \"nwa5\", \"nwa6\", \"cubanpete\", \"jimcarreymb\")]\n",
    "test_many = Tester(paths)\n",
    "test_many.use_negentropy()\n",
    "#test_many.use_kurtosis() # commented by default\n",
    "\n",
    "test_many.play_observed()\n",
    "test_many.mix_observed()\n",
    "test_many.play_observed()\n",
    "test_many.play_separated()\n",
    "#test_many.play_separated(sample_percent_of_data=0.5) # commented by default"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Konklusjon\n",
    "\n",
    "Vi har i denne notatboken implementert FastICA-algoritmen for å løse cocktailparty-varianten av blind kildeseparasjonsproblemet med uavhengig komponentanalyse.\n",
    "Separasjonsalgoritmen ble deretter testet på noen lydfiler av ulike antall, mengder og lengder ved bruk av både kurtose og negentropi.\n",
    "\n",
    "I nesten alle tilfeller ble lydfilene separert bemerkelsesverdig godt.\n",
    "Kun i ett tilfelle, der mange komponenter ble forsøkt separert med et tilfeldig utvalg datapunkter eller ved bruk av kurtose, feilet separasjonen ved at flere av sporene overlappet.\n",
    "Men til og med her gikk separasjonen feilfritt ved bruk av negentropi på alle datapunktene på maskinen vår.\n",
    "Kurtose og negentropi ser dermed ut til å fungere like godt og raskt til måling av gaussiskhet i de fleste tilfeller, men det virker som negentropi har muligheten til å oppnå større presisjon og kortere tidsbruk i enkelte mer krevende situasjoner.\n",
    "\n",
    "Vi har merket at FastICA-algoritmen kan benyttes svært effektivt til separasjon av få lange lydklipp når kun et tilfeldig utvalg av datapunktene benyttes i beregningen av demiksematrisen.\n",
    "Med et utvalg på kun 2% av datapunktene, separerte vi på denne måten fire 90-sekunders lydklipp feilfritt på under ett sekund på vår maskin!\n",
    "Legg merke til at den desidert mest tidkrevende operasjonen i denne testen var innlastingen av lydavspillingsverktøyene.\n",
    "Separasjon av de samme lydklippene med bruk av alle datapunktene produserte identiske resultater med betydelig større tidsbruk.\n",
    "Denne metoden fungerte derimot ikke like godt når vi forsøkte å benytte den på mange korte lydklipp.\n",
    "Her måtte vi benytte alle datapunktene for å få helt nøyaktige resultater, noe som tok betydelig lenger tid, selv om den totale lengden av lydklippene her var mye kortere.\n",
    "\n",
    "Vi har også sett at Python i samspill med biblioteket `numpy` er svært velegnet for å implementere en enkel og lærerik versjon av algoritmen.\n",
    "De effektive numeriske operasjonene i `numpy` og den enkle, tilgivende og intuitive syntaksen lot oss holde alt fokus på simpelthen å forflytte en matematisk beskrivelse av algoritmen over på datamaskinen.\n",
    "\n",
    "For å gjøre en bedre vurdering av FastICA-algoritmen, kunne vi testet den på et enda mer variert utvalg lydklipp av flere ulike antall, lengder og typer.\n",
    "Dette kunne avkreftet noen av konklusjonene vi har kommet frem til her, som kan være spesielle for de kombinasjonene av lydklipp vi har testet.\n",
    "Vi kunne også latt algoritmen gjette en tilfeldig initiell verdi for demiksematrisen, fremfor å alltid benytte den samme identitetsmatrisen, for å se om dette påvirker konvergensen.\n",
    "Vi kunne også forsket litt på å variere konvergenstoleransen i FastICA-algoritmen, men vi så ingen grunn til å gjøre dette her, da vi benyttet en svært lav toleranse i utgangspunktet og algoritmen alltid terminerte raskt nok når den ga korrekte separasjoner.\n",
    "\n",
    "Til tross for noen små imperfeksjoner når FastICA-algoritmen kjøres med noen spesielle innstillinger på visse kombinasjoner av lydklipp, kan vi uansett ikke gjøre annet enn å si at FastICA-algoritmen ser ut til å fungere utrolig godt til separasjon av mange, lange og varierte lydklipp!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Referanser\n",
    "[1] B. Owren. Blind source separation. 2019.\n",
    "\n",
    "[2] A. Hyvärinen, E. Oja. Independent component analysis: algorithms and applications. Neural Networks 12. 2000. Side 411-430."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
